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Abstract A general framework for the design of adaptive low-dissipative high order schemes 

is presented. It encompasses a rather complete treatment of the numerical ap- 
proach based on four integrated design criteria: (1) For stability considerations, 
condition the governing equations before the application of the appropriate nu- 
merical scheme whenever it is possible. (2) For consistency, compatible schemes 
that possess stability properties, including physical and numerical boundary con- 
dition treatments, similar to those of the discrete analogue of the continuum are 
preferred. (3) For the minimization of numerical dissipation contamination, ef- 
ficient and adaptive numerical dissipation control to further improve nonlinear 
stability and accuracy should be used. (4) For practical considerations, the nu- 
merical approach should be efficient and applicable to general geometries, and 
an efficient and reliable dynamic grid adaptation should be used if necessary. 
These design criteria are, in general, very useful to a wide spectrum of flow 
simulations. However, the demand on the overall numerical approach for non- 
linear stability and accuracy is much more stringent for long-time integration 
of complex multiscale viscous shock/shear/turbulence/acoustics interactions and 
numerical combustion. Robust classical numerical methods for less complex 
flow physics are not suitable or practical for such applications. The present ap- 
proach is designed expressly to address such flow problems, especially unsteady 
flows. The minimization of employing very fine grids to overcome the production 


* A chapter for Turbulent Flow Computation, (eds. D. Drikakis & B. Geurts), Kluwer Academic Publisher, 
2001. Expanded version of the paper for the Proceedings of the Third AFOSR International Conference on 
DNS/LES, Arlington Texas, August 4-9, 2001. Part of this work was carried out while the second author 
was a visiting scientist with RIACS, NASA Ames Research Center. 


1 


2 


of spurious numerical solutions and/or instability due to under-resolved grids is 
also sought [79, 17]. The incremental studies to illustrate the performance of 
the approach are summarized. Extensive testing and full implementation of the 
approach is forthcoming. The results shown so far are very encouraging. 

Keywords: Low-Dissipative Schemes, Adaptive Numerical Dissipation/Filer Controls, High 

Order Finite Difference Methods, Linear and Nonlinear Instabilities, Skew-Symmetric 
Form, Entropy Splitting, Summation-by-Parts, Integration-by-Parts, Wavelets, 
Multi-Resolution Wavelets, linear and nonlinear filters. 

1. Introduction 

Classical stability and convergence theory are based on linear and local lin- 
earized analysis as the time steps and grid spacings approach zero. This theory 
offers no guarantee for nonlinear stability and convergence to the correct solu- 
tion of the nonlinear governing equations. The use of numerical dissipation has 
been the key mechanism in combating numerical instabilities. Recent nonlinear 
stability analysis based on energy norm estimate [24] offers stability improve- 
ment in combating nonlinear instability for long-time integrations and/or simple 
smooth flows. These recent developments offer two major sources of stabilizing 
mechanisms, namely, from the governing equation level and from the numeri- 
cal scheme level. Employing a nonlinear stable form of the governing equation 
in conjunction with appropriate nonlinear schemes for initial-boundary-value 
problems (IB VPs) is one way of minimizing the use of numerical dissipations 
[56, 55]. On the other hand, even with the recent development, when em- 
ploying finite time steps and finite grid spacings in the long-time integration of 
multiscale complex nonlinear fluid flows, nonlinear instability, although greatly 
reduced, still occurs and the use of a numerical-dissipation/filter is unavoidable. 

Aside from acting as a post-processor step, most filters serve as some form of 
numerical dissipation. Without loss of generality, “numerical-dissipation/filter” 
is, hereafter, referred to as “numerical dissipation” unless otherwise stated. 
Proper control of the numerical dissipation to accurately resolve all relevant 
multiscales of complex flow problems while still maintaining nonlinear stabil- 
ity and efficiency for long-time numerical integrations poses a great challenge 
to the design of numerical methods. The required type and amount of numer- 
ical dissipation are not only physical problem dependent, but also vary from 
one flow region to another. This is particularly true for unsteady high-speed 
shock/shear/boundary-layer/turbulence/acoustics interactions and/or combus- 
tion problems, since the dynamics of the nonlinear effect of these flows are not 
well-understood [79], while long-time integrations of these flows have already 
stretched the limit of the currently available supercomputers and the existing 
numerical methods. It is of paramount importance to have proper control of 
the type and amount of numerical dissipation in regions where it is needed 
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but nowhere else. Inappropriate type and/or amount can be detrimental to the 
integrity of the computed solution even with extensive grid refinement. 

The present work is a sequel to [82, 83, 63, 65]. It is an expanded version of 
[84]. The objective here is to propose a rather complete treatment of the numer- 
ical approach based on the four integrated design criteria (i)-(4) stated in the 
abstract. The key emphasis here is to describe and illustrate with examples on 
an adaptive procedure employing appropriate sensors to switch on the desired 
numerical dissipation where needed, and leave the rest of the region free of 
numerical dissipation contamination, while at the same time improving nonlin- 
ear stability of the entire numerical process for long-time numerical integration 
of the complex multiscale problems in question. These sensors are capable of 
distinguishing shocks/shears from turbulent fluctuations and/or spurious high- 
frequency oscillations for a full spectrum of flow speeds and Reynolds numbers. 
The minimization of employing very fine grids to overcome the production of 
spurious numerical solutions and/or instability due to under-resolved grids is 
sought [17]. It was shown in [56, 20, 83, 63] that conditioning the govern- 
ing equations via the so called entropy splitting of the inviscid flux derivatives 
[83] can improve the over all stability of the numerical approach for smooth 
flows. Therefore, the same shock/shear detector that is designed to switch on 
the shock/shear numerical dissipation can be used to switch off the entropy 
splitting form of the inviscid flux derivative in the vicinity the discontinuous 
regions to further improve nonlinear stability and minimize the use of numerical 
dissipation. The rest of the sensors, in conjunction with the local flow speed 
and Reynolds number, can also be used to adaptively determine the appropriate 
entropy splitting parameter for each flow type/region. These sensors are readily 
available as an improvement over existing grid adaptation indicators [20], If 
applied correctly, the proposed adaptive numerical dissipation control is scheme 
independent, and can be a stand alone option for many of the favorite schemes 
used in the literature. 

Outline: A brief summary of linear and nonlinear stability and the logistics 
of advocating design criteria (l)-(4) for a complete numerical approach are 
discussed in Sections 2 — 4. Adaptive numerical dissipation controls for high 
order schemes are discussed in Section 5. Some representative examples to 
illustrate the performance of the approach are given in Section 6. 

2. Conditioning of the Governing Equations 

Traditionally, conditioning the governing partial differential equations (PDEs) 
usually referred to rewriting the governing equations in an equivalent set of 
PDEs in order to prove the stability and/or well-posedness of the PDEs. When 
numerical methods are used to solve PDEs that are nonlinear, it is well-known 
that different equivalent forms of the governing equations might exhibit dif- 
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ferent numerical stability, accuracy and/or spurious computed solutions, even 
for problems containing no shock/shear discontinuities. There are many con- 
ditioned forms of the governing equations proposed in the literature. Different 
conditioned forms of the nonlinear convection fluxes and the viscous fluxes 
have been proposed for the incompressible and compressible Navier-Stokes 
equations. Here we concentrate on the convection terms of these equations 
and mention a few conditioning forms which are precursors of the so called 
entropy splitting of the compressible Euler equations [83]. If a method-of-lines 
approach is used to discretize these equations, the entropy splitting reduces to 
the splitting of the convection flux derivatives. For the viscous terms, we only 
adapt the method of preventing odd and even decoupling on all of our numerical 
experiments whenever it is applicable [61]. 

The splitting of the nonlinear convection terms (for both the compressible 
and incompressible Navier-Stokes equations) into a conservative part and a 
non-conservative part has been known for a long time. In the DNS, LES 
and atmospheric science simulation literature, it is referred to as the skew- 
symmetric form of the momentum equations [4, 30, 5, 86]. It consists of the 
mean average of the conservative and non-conservative (convective form [86]) 
part of the momentum equations. The spatial difference operator is then ap- 
plied to the split form. From the numerical analysis standpoint, the Hirt and 
Zalesak’s ZIP scheme [27, 85] is equivalent to applying central schemes to 
the non-conservative momentum equations (convective form of the momentum 
equations). MacCormack [39] proposed the use of the skew-symmetric form 
for problems other than DNS and LES. Harten [25] and Tadmor [73] discussed 
the symmetric form of the Euler equations and skew-adjoint form of hyperbolic 
conservation laws, respectively. Although the derivation in these works is dif- 
ferent, the ultimate goal of using the split form is almost identical. This goal is 
to improve nonlinear stability, minimize spurious high-frequency oscillations, 
and enhance robustness of the numerical computations. The canonical splitting 
used by Olsson & Oliger [56] is a mathematical tool to prove the existence of a 
generalized energy estimate for a symmetrizable system of conservation laws. 
For the thermally perfect gas compressible Euler equations, the transformation 
consists of a convex entropy function that satisfies a mathematical entropy con- 
dition. The mathematical entropy function, in this case, can be a function of the 
physical entropy. Therefore, the resulting splitting was referred to as entropy 
splitting by Yee et al. [83]. The entropy splitting can be viewed as a more gen- 
eral form than its precursors which makes possible the L 2 stability proof of the 
nonlinear Euler equations with physical boundary conditions (BCs) included. 
The following subsections which were part of [65], provide more details. 
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2.1 Introduction to skew-symmetric splitting 

Consider a variable coefficient linear hyperbolic system 

Ut + A(x)U x = 0 a < x < b, (2.1) 

where U is a vector and the matrix A(x) is symmetric. Define the scalar product 
and norm, 

(U, V) = f U(xfV(x) dx, \\U\\ 2 = ( U , U). (2.2) 

J a 

It is possible to obtain an energy estimate for the solution by integration-by- 
parts. To do this, write the system in skew-symmetric form, 

U t + \{A(x)U) x + l -A{x)U x - l -A x U = 0. (2.3) 

Start from 

\f t W U W 2 = (^t) = -\[(U,(A(x)U) x ) 4- (U,A(x)U x ) + (UA X U)} 

(2.4) 

and perform the integration-by-parts 

([/, [A{x)U) x ) = ~(U x ,A(x)U) + [U T AU} b a = ~(U,AU X ) + [ U T AU ] b a , 

(2-5) 

where the last equality follows from the symmetry of A. This gives the energy 
norm estimate 

\j t W u W 2 = -\([U T AU} b a - (A X U,U)), (2.6) 

which is a standard result that has been known for a long time. It can be found 
in many textbooks on PDEs, e.g., [19]. 

For semi-discrete difference approximations, the same idea can be used. 
Introduce the grid points xj = a+(j — \)h,j = 1,2,... , TV on the interval [a, b\, 
with uniform grid spacing h = (b — a)/(N — 1). Apply a spatial discretization 
to the skew-symmetric form 

= A A {xj)DUj - l -D{A{xj)Uj) + \DA{xi)Uj (2.7) 

where D is a finite difference operator, approximating the spatial derivative. 
We will obtain an estimate in a discrete scalar product, 

N 

(I U,V) h = h'E (TijUiVj, 

ij—l 


( 2 . 8 ) 
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in the same way as for the continuous case. Here a hJ a positive definite matrix 
(identity matrix for the L 2 norm). The estimate becomes 

\i\\U\\l = -\{{U,ADU) h + (i U,D(AU )) h ) + \{U,D{A)U) h = 
-\{UlA N U N - U{ A X U\) + \{U,D{A)U) h , 

(2.9) 

where we now assume that the difference operator has the summation-by-parts 
(SBP) property, 

(U,DV) h = —{DU, V) h + U N V N - U x Vx (2.10) 

with respect to the discrete scalar product. The SBP property here is the discrete 
analogue of the integration-by-parts energy norm property. One simple SBP 
operator is given by 


DU\ = D-yU\ 

DUj = D 0 U V j = 2, 3, . . . , iV — 1 (2.11) 

DU n = D_U n 


for the scalar product 

h N ~ l h 

(i U,V) h = -UxVx +hJ2 UjVj + -U N V N . (2.12) 

^ i = 2 Z 

Here we define Dq Uj = (Uj+i — Uj-\)/{2h), D + = ( Uj + \ — Uj)/h,D^Uj = 

( Uj — Uj-\)/h. Higher order accurate SBP operators can be found; see [70]. 
For periodic problems, the SBP property is usually easy to verify. In this case 
the boundary terms disappear. 

The crucial point is the splitting of the convective term, 

AU X = l -AU x + 1 -{AU) X - l -A x U , (2.13) 

into one conservative and one non-conservative parts. The difference approxi- 
mation is applied to the split form. The skew-symmetric splitting for difference 
approximations has also been known for a long time. It was used in [30], and 
[31] to prove estimates for the Fourier method. See also [47], where SBP is 
proved for the Fourier method and a fourth-order difference method, when the 
boundaries are periodic. 

Although this L 2 estimate does not give uniform boundedness of the solu- 
tion, it has turned out in practical computations that methods based on skew- 
symmetric splitting perform much better for long-time integrations than un-spl it 
methods. 
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Actually, for a symmetric hyperbolic nonlinear system with A(x) replaced by 
A(U(x)), similar skew-symmetric splitting and energy norm can be obtained. 
One of the earliest works on skew-symmetric splitting is Arakawa [4], where a 
splitting was derived for the 2D Euler equations for incompressible fluid flow 
in vorticity stream function formulation. 


UJt tltyiOx (2.14) 

Here u> is the vorticity and ip is the stream function, such that the velocity (u, v) 
is (-ip y ,ip x ). In [4] it is shown that the approximation 

— — 5 (Dy , ip i Dx1>ijD y u>ij)+ 

\{D x {^i,jDyil>ij) — Dy(u>i t jD x ipi t j))-\- (2-15) 

\(Dy{lpi t jDx^i,j) — D x ('lp l j DytO l j ) ) 


leads to the estimates 


|IM| 2 = (^)/> = o 
{tp,ut t )h = 0. 


(2.16) 


Here it is assumed that boundary terms are equal to zero (homogeneous). D x 
and D v denote finite difference operators acting in the x- and j/ -direction re- 
spectively. The second estimate is related to the conservation of kinetic energy, 

^(IWf + IWf) = (217) 

The proof of the estimates only involves pairwise cancellation of terms accord- 
ing to the rule, 

(u, D(uv))h + (u,vDu)h = 0, (2.18) 

which holds for zero boundary data, if D satisfies (2.10). In [4], the operator 
(2.11) is used. Note that the use of u and ip in this section pertains to the 
vorticicity formulation symbols. In later sections, the same symbols will have 
different meanings. 

In [30], the inviscid Burgers’ equation, 

t it + (u 2 / 2) x = 0 (2.19) 

with the quadratic flux derivative split into \uu x + |('u 2 /2)x, is approximated 

= ~\( UjDuj + Du ))- 


as 


( 2 . 20 ) 
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For this approximation, we obtain 
\ d, 1 

2 ^IMI/i = ~^{(u,uDu) h + (■ u,Du 2 ) h ) = 0 (2.21) 

by using SBP on the last term. Again, boundary data is assumed to be zero. 
The split form (Du 2 + uDu )/ 3 is also used in [58]. 

It is important to note that the split form is non-dissipative in the sense that the 
highest grid frequency, Uj = ( — 1) J , has derivative zero, and thus cannot be 
seen or smoothed by the time integration. In addition, the difference operator 
applied to the split form should be done with D on the non-conservative term, 
and — D T , the negative adjoint of D, on the conservative term. However, since 
even order centered difference operators are anti-symmetric, we write D in 
both places. For odd orders of accuracy, the approximation should be done 
analogously to the first order example, 

^ = -I( UjD+Uj + D-uj). (2.22) 


2.2 Skew-Symmetric Splitting for Incompressible Fluid 
Flow 


For a 3-D incompressible Navier-Stokes equations of the form, 


u* + (u T V)u — — Vp + vAu, 
div u = 0, 


(2.23) 


skew-symmetric splitting can be applied on the convective terms to estimate 
the kinetic energy, u T u. Here the velocity vector is u = (ui,U 2 ^us) t , the 
pressure p and the viscosity coefficient is v. In [86] and in [26] the three forms 

(u r V)u (convective) 

^((u T V)u + div( uu T )) (skew-symmetric) (2.24) 

^(V(u r u) + u x V x u) (rotational) 

for the nonlinear terms are studied. They are equivalent to each other before 
the application of the numerical methods. Although the skew-symmetric and 
rotational forms are not in conservative form, they lead to conservation of kinetic 
energy that is important for long-time integration. For the inviscid case v — 0, 
when using the skew-symmetric form, we can estimate the kinetic energy as, 

^(IIV^II 2 ) = («1, («1 )t) + («2, («2)t) + (U3, (u 3 )t) = 

-( Ui , UjdjUi) - ( Ui , dj(uiUj)) - ( Ui , dip) = (2.25) 

~(ui,UjdjUi ) + ( djUi,UiUj ) + ( dtUi,p ) = 0 , 
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where the summation convention is used. We assume that boundary veloci- 
ties are zero. Integration-by-parts used above shows that the three pressure 
components are equal to (dzu(u),p), which is zero since the divergence is zero. 

In order to carry out the same estimate for a difference approximation, we 
should use the skew-symmetric form for the convective terms, and discretize 
by a difference operator having the SBP property. The convective terms then 
disappear from the estimate directly, without use of the divergence condition, 
just as for the PDE. To eliminate the pressure term, it is enough that the pressure 
derivatives and the divergence condition are discretized by the same difference 
operator (or by operators that are negative adjoints, in the case of odd order of 
accuracy). We end up with 


\f t \\ u L,k u i,j,k\\l < 0 (2.26) 

for the difference approximation in a discrete norm. Note that the inequality 
(2.26) should really be an equality in the setting it is proved. Perhaps it is 
possible to keep the viscous terms, and make it an inequality. 

The discrete estimate can also be derived from the discretized rotational 
form (as in [47]), but then the discretized divergence condition must be used 
to eliminate certain convective terms. For this to be possible, it is necessary 
that the divergence condition is discretized by the same operator as used for 
the other convection terms. Results from using the skew-symmetric form are 
compared with results from the rotational form in a turbulence simulation in 
[86]. The skew-symmetric form is found to give more accurate results. It is 
recommended that the pressure equation and the divergence condition should 
be discretized by the same SBP satisfying operator (or SBP operator for ease 
of reference), so that we can eliminate the term (u\,p x ) 4 - (u2,p y ) + (0,3, p z ). 

2.3 Skew-Symmetric Splitting for Compressible Fluid 
Flows 

Consider the equations of inviscid compressible fluid flow in one space di- 
mension 


[pn \ + ( pu^ + p\ = [o], (2.27) 

\ e ) t \u(e+p)J x \0 / 

with p, u , e and p, the density, velocity, total energy per unit volume and 
pressure, respectively. In [5], a skew-symmetric splitting of the convective 
terms in momentum equation is used. This splitting was originally presented 
in [16]. The discretization in [5] is made for a more general equation, but for 
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the simple equation (2.27), it becomes 

^ / Pj(t) \ / D(pjiij) 

— j ( pu)j(t ) 1 + I 2 D{PjUj) + 2 Pj U jDuj + iUjD(pjUj) + Dpj 

\ e?(t) / \ +Pj)) 

(2.28) 

where we now only consider semi-discrete approximations. The skew-symmetric 
splitting of the convective term in the momentum equation makes it possible to 
estimate 

d N N 

— ^ = — h ^2 u jDPj + boundary terms, (2.29) 

dt . 1 , 

3 = 1 J=1 

which is the discrete analogue of the estimate, 

— \\pu 2 \\ 2 = — (u,p x ) + boundary terms, (2.30) 

dt 

obtained from the PDE. 

2.4 Canonical (Entropy) Splitting for Systems of 
Conservation Laws 

The skew-symmetric splitting for the nonlinear incompressible and com- 
pressible Navier-Stokes equations discussed above only involve the nonlinear 
convective terms of the momentum equation, and not the entire inviscid flux 
derivatives of the PDEs. Actually, for a general nonlinear system of conserva- 
tion laws, 

U t + F{U) x = 0, a<x<b , 0<t, (2.31) 

we can perform a skew -symmetric splitting of the entire inviscid flux vector 
derivative F(U) X , if 

1 A(U) = dF/dU is symmetric. 

2 F is homogeneous, i.e., F(XU) = X^F(U), with (3 ^ —1. 

It is possible to show that 



A{U)U = 0F{U), ( 2 . 32 ) 

by differentiating the homogeneity relation with respect to A, and setting A = 1. 
Define the splitting as 


Ut + Y^A(U)U X + 


p 

1+/3 


F(U) X = 0 . 


( 2 . 33 ) 
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We then can show that 

\\U \\ 2 = (U,Ut) = -J^(U,A(U)U X ) - £_{U,F{U) x ) = 

A(U)U X ) + J X ,F(U)) - J ^{U T F{U)t = (2.34) 

-^r 0 [U T AU)l 

The two inner products disappear due to (2.32), and we are left with an estimate 
of the norm of the solution in terms of the solution on the boundary, a so 
called generalized energy estimate [56]. Olsson & Oliger [56] used the term 
generalized energy estimate because the norm will, in general, depend on an 
entropy vector W, but not the gradient of U in symmetrizable conservation 
laws. 

In [56] the splitting is first done without the homogeneity assumptions. In 
that case, the entropy flux function F E (U) 


F e (U)= C F(dU)d6 
Jo 

(2.35) 

is introduced. It follows that 


F/?(U)U = -F e + F, 

(2.36) 

so that 


F x = (F E U) X + F e U x . 

(2.37) 


If Fjj is symmetric, then F E is symmetric too, and we obtain, 


(U,F X ) = (U,(F E U) X ) + (U,F E U X ) = 

(^Z. jo ) 

~(U X ,F$U) + (U,F*U X ) + [UF$U} b a = [U T F§U] b a . 

If the flux function F is homogeneous, F E is just a scalar times F, and we 
recover the splitting (2.33). 

In practice, the Jacobian matrix A(U) is not symmetric, especially for more 
than 1-D. However, in many cases a symmetrizing variable transformation is 
available. The estimate for the homogeneous case above and symmetrizing 
transformations were given in [73]. In [56] the analysis is extended to non- 
homogeneous problems, and BCs are discussed in greater detail using the so 
called canonical splitting of the symmetrizable conservation laws. Formulas 
for symmetrizing the nonlinear compressible Euler equations are given in [25], 
and the corresponding analytical form for the canonical splitting of the perfect 
gas compressible case is given in [20]. It was further extended to thermally 
perfect gases and to 3-D generalized coordinates that preserve freestream in 
[83,77]. 
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Let the symmetrizing vector W be related to U via a transformation, 

U = U{W). (2.39) 

It can be proved that the existence of a symmetrizing transformation is equiva- 
lent to the existence of an entropy function, E(U ), with Euu positive definite. 
The entropy function is a function such that 

E t + F x e = 0 (2.40) 

is an additional conservation law, obtained by multiplying the original conser- 
vation law by Ejj(U). F E is the entropy flux, related to the entropy by 

(F e ) t = ElFu , (2.41) 

an equation which is overdetermined, and therefore does not have a solution 
for all systems of conservation laws. Written in terms of the entropy function, 
(the inverse of) the symmetrizing change of variables (2.39) is defined by W = 
Eu(U). The change of variables dU /dW is symmetric and positive definite, 
and the new Jacobian dF/dW is symmetric. If furthermore U and F are 
homogeneous in W of degree /?, which is the case for the thermally perfect gas 
Euler equations for any (3^—1, the formulas become simple. In that case we 
insert the change of variables into the conservation law and obtain 

U w W t + F w W x = 0, (2.42) 

and define the split form of the flux derivative [56] 

Ut + + TTp FwW * = °> (2 - 43) 

with (3 a splitting parameter (/? = oo recovers the original conservative form). 
Here (3^—1 and, for a perfect gas, (3 > 0 or f3 < y— ;. The theory only 
gives the range of f3 and does not give any guidelines on how to choose (3 for 
the particular flow. The vectors Fw and W can be cast as functions of the 
primitive variables (p, u,p) and (3. From the study of [83], (3 is highly problem 
dependent. Multiplying the above equation by W and integrating gives 

-(1 + P)(w, Ut) = P(w, F x ) + (W, F w W x ) = 0(W, F x ) + ( F w W , W x ). 

(2.44) 

Integration-by-parts in space gives 

(1 + 0)(W, U t ) = ~[W T F w W] b a . (2.45) 

We thus obtain the estimate 

&(W, UwW) = (Wt, UwW) + (W, (U w W)t) = (Ut, W) + (5(W, U t ) = 
(1 +P)(W,U t ) = -[W T F w W) b a . 


(2.46) 
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In order to have an energy estimate, the boundary term [ It ' / F\ t - W ] should be 
of the sign that makes the time derivative of the norm negative. For stability 
the entropy norm (W, UwW) should be bounded. 

It is noted that the energy estimate can be shown to be identical with 

| J E(U) = [F E ) b a (2.47) 

obtained by integrating the entropy equation (2.40) in space. It follows that, 
WUwW = (1 + (S)E(U) [73]. We can show that for the thermally perfect 
gas compressible Euler equations the mathematical entropy function can be a 
function of the physical entropy. Therefore, the resulting splitting was referred 
to as entropy splitting by Yee et al. [83]. 

3. Discrete Analogue of the Continuum 

Standard stability guidelines for finite difference methods in solving non- 
linear fluid flow equations are based on a linearized stability analysis. The 
linear stability criterion is applied to the frozen nonlinear problem at each time 
step and grid point. Most often the numerical BC (or boundary scheme), if 
needed, is not part of the stability analysis. The preceding section summa- 
rizes the historical perspective of entropy splitting of the fluid flow equation 
related to stable spatial finite discretization without paying attention to numer- 
ical BC. This section expands on stable finite difference methods that have a 
discrete analogue of the conditioned governing equations IB VPs. For ease of 
reference, “scheme” or more precisely “interior scheme” here generally refers 
to spatial difference schemes for the interior grid points of the computational 
domain, whereas “boundary scheme” is the numerical boundary difference op- 
erator for grid points near the boundaries. However, without loss of generality, 
we also adopt the conventional terminology of denoting “scheme” interchange- 
ably as either the “combined interior and boundary scheme” or just the “interior 
scheme” within the context of the discussion. 

The only tool needed to derive the norm estimates presented in the preceding 
section was integration-by-parts. One main point in this section is that the same 
norm estimates can be made for a semi-discrete difference approximations, if 
the differential operators are approximated by difference operators having the 
SBP property. Examples of this were shown in (2. 16) and (2.21). Of course, the 
other estimates in the preceding section can be carried over to a semi-discrete 
approximation by use of SBP difference operators. 

The discussion is divided into linearly stable and nonlinearly stable differ- 
ence methods. It is important to point out that when solving the Navier-Stokes 
equations with complex viscous shock, shear-layer, and boundary layer and/or 
chemical reaction interactions, even after incorporating tools from recent de- 
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velopments, the finite difference methods considered, although more rigorous 
than standard algorithms, are only linearly stable in a strict sense. 

3.1 Linearly Stable Difference Methods 

The most basic linear stability criterion is to investigate the behavior of the 
difference method when applied to a problem of constant coefficients and peri- 
odic boundaries. The Fourier symbol of the operator should be bounded. For 
higher than first-order methods, a complication is introduced by the numerical 
boundary treatment. Norm estimates, or normal mode analysis are normally 
employed (See Gustafsson et al. [24]). With these methods it is possible to 
prove stability for linear IBVPs. Difference operators having the SBP prop- 
erty with numerical BCs included have recently received some attention. See 
Strand, Olsson, and Nordstrom & Carpenter [70, 53, 54, 55, 51]. As already 
discussed in Section 2, the idea with these operators is to have the property 

(. DU 3 , Vj) h = ~{Uji DVj) h + V N U N - V X U U (3.1) 

where D is a difference operator approximating d/dx , including the accom- 
panied boundary scheme. Typically D is a standard centered operator in the 
interior of the computational domain, and has a special one-sided form near 
boundaries. The discrete scalar product is defined by (2.12), and is weighted 
by a positive definite matrix, a. For the standard L 2 -norm, a is the identity 
matrix. In [53, 70], formulas for the norm and boundary modifications of D 
are given which ensure the SBP property for operators up to order of accuracy 
eight. SBP satisfying numerical BCs are very different from the traditional 
numerical boundary treatment. For example, for a sixth-order central interior 
scheme, the SBP satisfying boundary schemes involved the modification of the 
central scheme at least 6 points from the boundary. The coefficients of these 
SBP boundary schemes are rational and irrational fractions. The coefficients 
of the boundary scheme are determined together with the weights, cr, in the 
scalar product, so that for each operator there is a particular scalar product in 
which the SBP property holds. With the SBP property, norm estimates of the 
difference approximation can be accomplished as the discrete analogue of the 
continuous energy estimate of the corresponding IBVP of the PDE. 

3.2 Nonlinearly Stable Difference Methods 

When using a linearly stable method on a nonlinear problem, nonlinear in- 
stabilities can appear. Instabilities can appear already for a linear problem with 
variable coefficients. For variable coefficient problems, it can be proved that 
numerical dissipation of not too high order will make the method stable. From 
a theorem by Strang it follows that a finite difference approximation of a non- 
linear problem is stable, if the variable coefficient linearized approximation is 
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stable, and the solution and the difference scheme are smooth functions. This 
is one reason for using numerical dissipation in practical flow simulation [24]. 

Before 1994, rigorous stability estimates for accurate and appropriate bound- 
ary schemes associated with fourth-order or higher spatial interior schemes were 
the major stumbling block in the theoretical development of combined interior 
and boundary schemes for nonlinear systems of conservation laws. Olsson 
[55] proved that an energy estimate can be established for second-order central 
schemes. To obtain a rigorous energy estimate for high order central schemes, 
one must apply the scheme to the split form of the inviscid governing equation. 
A discrete analogue of the continuum using a semi discrete approach can be 
written as 

f = <3 - 2) 

Here, D is a difference operator, having the SBP property [55, 70]. The estimate 

j t (W, UwW) h = -Wj,F w {W N )W N + WfFwiW^W! (3.3) 

in the discrete scalar product follows in the same way as for the PDE with 
indices 1 and N the end points of the computational domain, and h the grid 
spacing. Here the SBP satisfying difference operator, for example, consists of 
central difference interior operators of even order together with the correspond- 
ing numerical boundary operators that obey the discrete energy estimate. See 
Olsson and Strand for forms of the SBP boundary operators [55, 70]. As noted 
in Section 2, if odd order of the spatial discretizations are used, the difference 
operator D in (3.2) should be modified. In this case, D should be employed on 
the non-conservative term, and — D T , the negative adjoint of D on the conser- 
vative term, i.e., 

- -Th { - DT)F{U ’ } - TTp mVi)OWi - (3 ' 4) 

3.3 SBP Difference Operators and Full Discretization 

There are two additional difficulties when applying the above semi-discrete 
SBP spatial discretization methods to realistic problems. 

■ How to impose given physical BCs without destroying the SBP property. 
For example, assume that we are given a boundary value u(x i) = g 
at the leftmost grid point of the domain at j = 1. Applying the SBP 
operator at j = 2, 3, . . . , N, and imposing u\ = g would not lead to an 
estimate, since the one sided operator that should have been applied at 
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j — 1 is not applied. This in turn leads to additional non-zero terms in 
the scalar product (e.g., from the operator at j = 2) which should have 
been canceled by terms from the operator at j = 1. 

■ How to discretize in time. The semi-discrete estimates that show the 
norm decreases, will not necessarily lead to a decrease of the norm in 
a time discretized approximation. In many practical cases we do obtain 
a small increase in the norm for the time discrete problem. However, 
the stability is still greatly improved by use of entropy splitting and SBP 
operators, when compared with more standard schemes. 

Several methods showing how to impose the physical BC have been proposed 
to overcome the first difficulty. Examples are the projection method [53] and the 
penalty method called “simultaneous approximation term” (SAT) [10, 51]. For 
comparison of these methods, see [69, 48, 29]. The methods given in [10, 51] 
and [53] are based on linear properties and cannot be trivially generalized to 
the nonlinear Navier-Stokes equations, except for certain special cases. One 
such special case where the nonlinear case is covered by the theory involves 
imposing velocity zero on solid walls, where the simple approach of setting the 
velocity to zero after each time step coincides with the projection method in 
[53]. 

In addition, when time-dependent physical boundaries are involved, an ad- 
ditional complication arises, especially for multi-stage Runge-Kutta methods. 
If the time-dependent physical BC is not imposed correctly, the overall order 
of accuracy of the scheme cannot be maintained. Some systematic remedies 
are proposed but are rather complicated for variable coefficients and even more 
complicated for nonlinear problems. See [10, 8, 29]. 

For the full discretization of the problem, we should discretize in time in 
such a way that the discrete energy estimate also holds. The obvious solution 
would be to discretize in time in a skew-symmetric way, in a manner similar to 
the spatial discretization, e.g., 

ih D ‘ u i «7= (35) 

-£fDF(VJ) - ^F w (WV)DWf, 

where D t is a difference operator acting in the time direction. However, it turns 
out that the SBP property of the time difference quotient leads to a problem 
which is coupled implicitly in the time direction. To solve it we have to solve 
a nonlinear system of equations for all time levels in the same system, leading 
to an impractically large computational effort. Furthermore, numerical experi- 
ments shown in Sjogreen & Yee [65] indicated that a bounded L 2 entropy norm 
(VF, UwW)h does not necessarily guarantee a well behaved numerical solution 
for long-time integrations. In other words, L 2 stability does not necessarily 
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guarantee an accurate solution. In practical computations, the classical Runge- 
Kutta time discretizations using the method-of-lines approach (which we used 
for our numerical experiments) works well, but we have not been able to prove 
a time discrete entropy estimate for this method. In addition, numerical exper- 
iments shown in [65] indicate that the time discrete problem does not have a 
decreasing entropy norm for all values of (3. Numerical experiments in Yee et 
al. [83, 61] also indicate the wide variations of the (3 value for a full spectrum 
of flow problems. For example, if a constant j3 is used for problems containing 
shock waves, a very large value of / 3 is needed. Otherwise divergent solutions 
or wrong shock location and/or shock strengths are obtained. In view of these 
findings, employing a constant (3 (within the allowable range of (3) throughout 
the entire computational domain appears not to be the best approach unless the 
flow problem is a simple smooth flow. Studies in [83, 63] indicate that the split 
form of the inviscid flux derivatives does help in minimizing the use of numer- 
ical dissipation. What is needed is adaptive control of the (3 parameter from 
one flow region to another as well as from one physical problem to another. We 
caution that if the adaptation is not handled correctly, an abrupt switching of 
the (3 can introduce spurious jumps in the numerical solutions. See [65] for the 
discussion. 

In our computer code for the numerical experiment, we have implemented 
the sixth-order SBP operators by the projection and SAT methods given in [10, 
51, 53]. They both perform satisfactorily, and no big difference in performance 
has been observed between them. See Section 6 for a 3-D compressible channel 
flow computation. We note that the majority of the physical boundaries of our 
viscous models are not time dependent, and the loss of spatial accuracy due to 
the multistage Runge-Kutta method is not a major concern. 

4. Adaptive Numerical Dissipation Control 

This section discusses the need for adaptive numerical dissipation controls 
in addition to conditioning the governing equations. An advanced numerical 
dissipation model for multiscale complex viscous flows is described. 

The linear and nonlinear numerical dissipation (not filter) presently available 
is either built into the numerical scheme or added to the existing scheme. The 
built-in numerical dissipation schemes are, e.g., upwind, flux corrected transport 
(FCT), total variation diminishing (TVD), essentially non-oscillatory (ENO), 
weighted ENO (WENO), and hybrid schemes (e.g., those that switch between 
spectral and high-order shock-capturing schemes). The built-in nonlinear nu- 
merical dissipation in TVD, ENO and WENO schemes was designed to capture 
accurately discontinuities and high gradient flows while hoping to maintain the 
order of accuracy of the scheme away from discontinuities. These schemes have 
been shown to work well in a variety of rapidly developing shock-shock inter- 
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actions that do not involve multiscale physics or long-time wave propagations. 
For multiscale physics that require low-dispersive errors, the amount of numer- 
ical dissipation built into these schemes is not optimal. In addition, analogous 
SBP theory for these schemes is not available. Moreover, they are more com- 
putationally expensive than standard high order centered schemes, and have 
severe limitations on the order of accuracy in the vicinity of the discontinu- 
ities and steep gradient regions. The inaccuracy of the numerical solutions can 
contaminate the entire flow field downstream [14]. Although, the amount of 
numerical dissipation is less than linear numerical dissipations, when applied 
to convection portions of viscous flows, it conflicts with the physical viscosity 
and can wash out the true physical steep gradient and/or turbulent structures. 
Aside from this fact, viscous reacting flows are even more difficult to simulate 
than non-reacting viscous flows. In the presence of numerical dissipations, even 
what is believed to be the optimal amount for non-reacting flows might have 
detrimental effects, e.g., wrong speeds of propagation and/or spurious traveling 
waves [36, 33, 34]. 

There exist different specialized linear and nonlinear filters to post process 
the numerical solution after the completion of a full time step of the numerical 
integration. Since they are post processors, the physical viscosity, if it exists, is 
taken into consideration. The main design principle of linear filters is to improve 
nonlinear stability, to stabilize under-resolved grids [17] and to de-alias smooth 
flows, while the design principle of nonlinear filters is to improve nonlinear 
stability as well as accuracy near discontinuities. When discontinuities are 
present in the solution, linear filtering and/or entropy splitting might not be 
helpful or not applicable. The nonconservative terms of the entropy splitting 
might lead to inconsistent behavior at shocks/shears [83, 63, 65]. See, for 
example, [22, 17, 18, 76] for forms of linear filters, and see [82, 83, 63] for 
forms of nonlinear filters. The use of the linear filter concept for smooth and/or 
turbulent flows has been employed for over two decades [76, 2, 37, 18]. For 
direct numerical simulation (DNS) and large eddy simulation (LES), there are 
additional variants of the linear filter approach. It was shown in Fischer & 
Mullen [17] that adding an appropriate filter can stabilize the Galerkin-based 
spectral element method in convection-dominated problems. The Fischer & 
Mullen numerical example illustrates the added benefit of the high-order linear 
filter. See Section 5.5 for a discussion. 

For the last decade, CPU intensive high order schemes with built-in non- 
linear dissipation have been gaining in popularity in DNS and LES for long- 
time integration of shock-turbulence interactions. Aside from the aforemen- 
tioned short-coming of these built-in nonlinear dissipation high order schemes, 
their flow sensing mechanism is not sophisticated enough to clearly distinguish 
shocks/shears from turbulent fluctuations and/or spurious high-frequency oscil- 
lations. In [82, 83, 63] it was shown that these built-in numerical dissipations 
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are more dissipative and less accurate than the nonlinear filter approach of 
[82, 83, 63] with a similar order of accuracy. It was also shown that these 
nonlinear filters can also suppress spurious high-frequency oscillations. How- 
ever, a subsequent study of Sjogreen & Yee [65] showed that the high order 
linear filter can sustain longer time integration more accurately than the non- 
linear filter for low speed smooth flows. In other words, for the numerical 
examples that were examined in [65], the high order linear filter can remove 
spurious high-frequency oscillation producing nonlinear instability better than 
the second-order nonlinear filter. Higher than third-order nonlinear filters might 
be able to improve their performance or might outperform the high order lin- 
ear filters in combating spurious high-frequency oscillations at the expense 
of more CPU time and added complexity near the computational boundaries. 
These findings prompted the design of switching on and off or blending of 
different filters to obtain the optimal accuracy of high order spatial difference 
operators as proposed in Yee et al. and Sjogreen & Yee [83, 63]. The missing 
link of what was proposed in [83, 63] is an efficient, automated and reliable 
set of appropriate sensors that are capable of distinguishing key features of the 
flow for a full spectrum of flow speeds and Reynolds numbers. 

We propose to enhance the conditioning of the equations with an advanced 
numerical dissipation model, which includes nonlinear sensors to detect shocks/ shears 
and other small scale features, and spurious oscillation instability due to under- 
resolved grids. Furthermore, we will use the detector to switch off the entropy 
splitting at shocks/shears and adjust the entropy splitting parameter with the aid 
of the local Mach number and Reynolds number in smooth regions as discussed 
earlier. The advanced numerical dissipation model can be used: (Option I) as 
part of the scheme, (Option II) as an adaptive filter control after the completion 
of a full time step of the numerical integration or (Option III) as a combination 
of Options I and II. For example, we can combine high order nonlinear dissi- 
pation (with sensor control) using Option I and nonlinear filter (with a different 
sensor control) using Option II. 

The numerical experiments we have conducted so far concentrate on an 
adaptive procedure that can distinguish three major computed flow features to 
signal the correct type and amount of numerical dissipation needed in addi- 
tion to controlling the entropy splitting parameter. The major flow features 
and numerical instability are (a) shocks/shears , (b) turbulent fluctuations, and 
(c) spurious high-frequency oscillations. It is important to not damp out the 
turbulent fluctuations. The procedure can be extended if additional refinement 
or classification of flow types and the required type of numerical dissipation 
is needed. There exist different detection mechanisms in the literature for the 
above three features. These detectors are not mutually exclusive and/or are 
too expensive for practical applications. We believe that the multiresolution 
wavelet approach proposed in Sjogreen & Yee [63] is capable of detecting all 
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of these flow features, resulting in three distinct sensors. If chosen properly, 
one multiresolution wavelet basis function might be able to detect all three fea- 
tures. For an optimum choice, one might have to use more than one type of 
wavelet basis functions but at the expense of an increase in CPU requirements. 
Some incremental studies into the use of entropy splitting and the application 
of these sensors were illustrated in [83, 63, 64, 66, 65, 50]. The next section 
summarizes the development of adaptive filters for a special class of high order 
discretizations. 

5. High Order Filter Finite Difference Methods 

The adaptive numerical dissipation controls discussed in the preceding sec- 
tion are scheme independent. This section applies the adaptive ideas to high or- 
der central schemes. We first summarize the high order nonlinear filter schemes 
that were developed for shock/shear capturing. We then extend these filtering 
ideas to include more complex flow structures by blending more than one filter 
and sensing tool. 

5.1 ACM and Wavelet Filter Schemes for Discontinuity 
Capturing 

An alternative to linear filter and/or standard shock-capturing schemes for 
viscous multiscale and long-time wave propagation computations is the ACM 
(artificial compression method) and wavelet filter schemes described in [82, 
63]. A high order centered base scheme together with the nonlinear dissipative 
portion of a shock-capturing scheme, activated by an ACM or wavelet sensor 
is used as the filter. Often an entropy split form of the inviscid flux derivatives 
is used. The idea of the ACM filter scheme is to have the spatially higher non- 
dissipative scheme activated at all times and to add the full strength, efficient 
and accurate numerical dissipation only at the shock layers and steep gradients. 
Thus, it is necessary to have good detectors which flag the layers, and not 
the oscillatory turbulent parts of the flow field. While minimizing the use 
of numerical dissipation away from discontinuities and steep gradients, the 
ACM filter scheme consists of tuning parameters that are physical problem 
dependent. To minimize the tuning of parameters, new sensors with improved 
detection properties were proposed in Sjogreen & Yee [63]. The new sensors 
are derived from utilizing appropriate non-orthogonal wavelet basis functions, 
and they can be used to completely switch off the extra numerical dissipation 
outside shock layers. The non-dissipative spatial base scheme of arbitrarily high 
order of accuracy can be maintained without compromising its stability at all 
parts of the domain where the solution is smooth. The corresponding scheme 
is referred to as the wavelet filter scheme. This nonlinear filter approach is 
particularly important for multiscale viscous flows. The procedure takes the 
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physical viscosity and the reacting terms into consideration since only non- 
dissipative high order schemes are used as the base scheme. In other words, 
numerical dissipations based on the convection terms are used to filter the 
numerical solution at the completion of the full step of the time integration at 
regions where the physical viscosity is inadequate to stabilize the high frequency 
oscillations due to the non-dissipative nature of the base scheme. 

The method applied to the 2-D conservation law where U is the conservative 
vector and F and G are the inviscid fluxes, 

U t + F(U) x + G(U) y = 0, (5.1) 

can be described as taking, e.g., one full time step by a Runge-Kutta method on 
the semi discrete system without or with entropy splitting, respectively, by 

^ = —DjF(Uj t k) - D K G{U hk ), (5.2) 

dJ ¥ =-T U D jF{U h k) + D K G{U hk )\ (53) 

TT0[F W (U hk )DjW^ k + G w (U jjk )D K W Jtk }, 

where Dj and Dk are high order finite difference operators, acting in the j- and 
fc-direction respectively. They can be the SBP satisfying higher-order difference 
operators (e.g., sixth-order central scheme with SBP boundary schemes). We 
consider here a rectangular grid with grid spacing Ax and Ay and time step 
At. Denote a full Runge-Kutta step by 

Uji 1 = RK(U^ k ). (5.4) 

After the completion of a full Runge-Kutta step, a filter (post processing) step 
is applied leading to 

Ufi 1 = Uji 1 - A x(F j+l/2ik - Fi_ w ) - A y(G jM1/2 - G Jtk _ 1/2 ) (5.5) 

with A x = At/ Ax and \ y = At/ Ay. The filter numerical fluxes Fj + i/ 2 ,k 
and Gj yk + 1/2 act in the j- and k- coordinate directions respectively, and are 
evaluated on the function U n+l . For example, 

Fj+ 1/2, k = 2' R J + 1 /2^f+l/2 ( 5 - 6 ) 

where Rj+ 1/2 is the right eigenvector matrix of the Jacobian of the inviscid flux 
F evaluated at Roe’s average state with the k index suppressed. The Zth element 
of the filter flux d>* +1 / 2 in the a; -direction (4> l j+i/ 2 )* is a product of a sensor 
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u>j +1 /2 an d a nonlinear dissipation < /A +1 j 2 , l — 1, 2, 3, 4. With the omission of 
the k index, it is of the form 

(4> l j+l/ 2 )* — a ’j+l/ 2 < ^j'+l/ 2 - (5-7) 

For the ACM sensor, wj +1/2 is a product of a physical dependent sensor coeffi- 
cient and a gradient like detector. The nonlinear numerical dissipation ^j+ 1/2 
can be obtained, from the dissipative portion of a TVD, ENO or WENO scheme. 
For example, the numerical flux i/ j+ i/ 2 ,fc of a second- or third-order TVD, ENO 
or WENO scheme can be written 

H j+ 1 / 2 , fc = \{F hk + F j+ i, fc ) + ^+i/2*j+i/2, (5-8) 

with the first two terms corresponding to the flux average of a centered difference 
and $j+i /2 with elements <t> l j+1 / 2 being the numerical dissipation portion of 
the scheme. 

For all the numerical experiments, the numerical dissipation portion of the 
Harten-Yee scheme is used. It has the form for the j-direction 


<£5+1/2 = ^(«5+1/2)(s5+1 + 9j) - Q( a j+l/2 + ^ + l/2) 5 5+l/2 (5 ' 9) 


with Q{x) = \/x 2 -I- e 2 , the entropy satisfying remedy for the scheme with 
entropy correction parameter e (not to be confused with the entropy splitting 
parameter). a* +1/ , 2 *s the Ith characteristic speed evaluated at the Roe’s average 

state in the j-direction. 7 y + i /2 is the modified characteristic speed and g l j is 
a slope limiter which is a function of the jump in the characteristic 

variable in the x -direction. 

A form of the ACM sensor w' +1/2 proposed in [82] is 

uj l j+ 1/2 = Kmax(9pQ l j + i) ( 5 . 10 ) 


where 


0j = 


l a 5+l/2l 155-1/2 
i 5 5+i /2 i + |5 5-i/2 


(5.11) 


See [82, 83] for details. It was shown in [63] that the method can be improved by 
letting the sensor w( +1/ , 2 be based instead on a regularity estimate obtained from 
the wavelet coefficients of the solution. The wavelet analysis gives an estimate 
of the so called local Lipschitz exponent a. The dissipation is switched on for 
low a values, and switched off when a becomes large [63]. The wavelet analysis 
is more general and can be used to detect other features besides shocks/shears. 
The following gives a more detailed explanation. 
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5.2 Wavelet Sensor for Multiscale Flow Physics 

Wavelets were originally developed for feature extraction in image process- 
ing and for data compression. It is well known that the regularity of a function 
can be determined from its wavelet coefficients [13, 46, 41] far better than from 
its Fourier coefficients. By computing wavelet coefficients (with a suitable set 
of wavelet basis functions), we obtain very precise information about the regu- 
larity of the function in question. This information is obtained just by analyzing 
a given grid function. No information about the particular problem which is 
solved is used. Thus, wavelet detectors are general, problem independent, and 
rest on a solid mathematical foundation. 


As of the 1990’s, wavelets have been a new class of basis functions that 
are finding use in analyzing and interpreting turbulence data from experiments. 
They also are used for analyzing the structure of turbulence from numerical data 
obtained from DNS or LES. See Farge [15] and Perrier et al. [57]. There are 
several ways to introduce wavelets. One standard way is through the continuous 
wavelet transform and another is through multiresolution analysis, hereafter, re- 
ferred to as wavelet based multiresolution analysis. Mallet and collaborators 
[41, 42, 43, 44, 45, 46] established important wavelet theory through multires- 
olution analysis. See references [72, 71] for an introduction to the concept of 
multiresolution analysis. Recently, wavelet based multiresolution analysis has 
been used for grid adaptation (Gerritsen & Olsson [20]), and to replace existing 
basis functions in constructing more accurate finite element methods. Here we 
utilize wavelet based multiresolution analysis to adaptively control the amount 
of numerical dissipation. 

The wavelet sensor estimates the Lipschitz exponent of a grid function fj 
(e.g., the density and pressure). The Lipschitz exponent at a point x is defined 
as the largest a satisfying 


\f(x + h) - f(x)\ 

So h- 


<C, 


(5.12) 


and this gives information about the regularity of the function /, where small 
a means poor regularity. For a C 1 wavelet function ip with compact support, 
a can be estimated from the wavelet coefficients, defined as 

WrnJ =< >= J f{x)tpm,j{x)dx, (5.13) 

where 

^• = 2 ">(^) (5 - 14) 
is the wavelet function xp m j on scale m located at the point j in space. This 
definition gives a so called redundant wavelet, which gives (under a few tech- 
nical assumptions on a non-orthogonal basis for L 2 . Theorem 9.2.2 in [13] 
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states that if \p is C 1 and has compact support, and if the wavelet coefficients 
in a neighborhood of jo decay as 2 ma as the scale is refined, then 


max ? \w. 


m,j 


the grid function fj has Lipschitz exponent a at jo- In practical computations, 
we have a smallest scale determined by the grid size. We evaluate w m j on this 
scale, mo, and a few coarser scales, mo + 1, mo 4- 2, and estimate the Lipschitz 
exponent at the point jo by a least square fit to the line [63] 


max logo I w m A = mot + c. 
j near jo 621 J 1 


(5.15) 


Proper selection of the wavelet ip is very important for an accurate detection 
of the flow features. The result in [46, 45], which is used in [20], gives the 
condition that ip(x) should be the /cth-derivative of a smooth function rj(x) 
with the property 


rj(x) > 0, 


J rj{x) dx 


1 


lim r? (fc) (a;)= 0 . 

x — »±oo 


(5.16) 


Then the result is valid for 0 < a < k. A continuous function f(x) has a 
Lipschitz exponent a > 0. A bounded discontinuity (shock) has a = 0, and a 
Dirac function (local oscillation) has a = — 1. Large values of A; can be used in 
turbulent flow so that large vortices or vortex sheets can be detected. Although 
the theorem above does not hold for a negative, a useful upper bound on a can 
be obtained from the wavelet coefficient estimate. 

For the numerical experiments, the wavelet coefficient w m j is computed 
numerically by a recursive procedure, which is a second-order B-spline wavelet 
or a modification of Harten’s multi-resolution scheme [63]. We can express the 
algorithm as follows. Introduce the grid operators 


Afj — 12k=~p dkfj+k 

Dfj ~ 5^/c-_ p c kfj+k 

and its mth level expanded versions 

Amfj — d>kfj+2 rn k 

D m fj = Y^k--p c kfj+2 m k-> 


(5.17) 


(5.18) 


where the integers p and q and the coefficients c & and are related to the 
chosen ip( x ) and </>(x), and can be determined from them. Here cj)(x ) is the so 
called scaling function of the multiresolution wavelets. 

The mth level of wavelet coefficients can be written as 


%,j — (/ ^ VVnj) — • • ■ ^0 fj 5 ^Tl — 1,2,.... (5.19) 
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Once the coefficients <4 and c/~ are determined, the computation is a very 
standard application of grid operators. In practice, we only use mo — 3 to 5. 
To be able to compute up to the boundary, we use one sided versions of the 
given operators. This seems to work well in practice, although it is not covered 
by the wavelet framework described above. 


Detectors from the B-Spline Wavelet Basis Function. Developing the 
best suited wavelets that can characterize all of the flow features might involve 
the switching or blending of more than one mother wavelet 'ip(x) and scaling 
function 0(x), especially if one needs to distinguish turbulent fluctuations from 
shock/shear and/or spurious high frequency oscillations. The mother wavelet 
function used in [20] and described in detail in [46] meets some of our require- 
ments. It is obtained from second order B-splines. 


f0 


'ip(x) = < 


—2(x — l) 2 
— 4x(l — x) + 2x 2 
— 4x(l + x) — 2x 2 
2(x + l) 2 
0 


x > 1 

1/2 < x < 1 
0 < x < 1/2 
— 1/2 < x < 0 
-1 < x < -1/2 
x < — 1 


For this wavelet (5.20), there exists a scaling function, given by 




'0 x > 2 

\{x — 2) 2 1 < x < 2 

< -(x- 1/2) 2 + 3/4 0 < x < 1 

^(x + l) 2 — 1 < x < 0 

i 0 x < — 1 


(5.20) 


(5.21) 


The normalization is such that the integral of the scaling function above is 
equal to one. The functions above are standard, and can be found in [13]. The 
scaling function differs by a shift from the scaling function used in [20], but the 
important relations 

4>{x) = \<f>{2x + 1) + f <M2x) + |0(2x - 1) + \4>{2x - 2) 

ip(x) = <j)(2x + 1) — 4>{ 2x) 

hold, and give the grid operators 

Afj = {fj-i + 3 fj + 3/j+i + fj + 2)/8, j = 2, ... ,N — 2 ^23) 

Dfj )/2 j = 2, . . . ,N. 

Note that this wavelet stencil is not symmetric. In general, the wavelet 
coefficients involve points from p2 mo ~ 1 to —q2 m °~ 1 , giving a stencil of totally 
(p + q) 2 mo_1 + 1 points. 
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Detectors from the Redundant Form of Harten Multiresolution Wavelet. 

For the redundant form of Harten multiresolution wavelet there is more than 
one choice for the interpolation function. See Sjogreen [62] for a discussion. 
The exact form of the method for the computations in this article is 

Af, =(/,--! + /;+ 1)/2 j = 2,...,N-l (524) 

Dfj = fj - Afj j = 2, . . . , N - 1. 

The above choice was made in order to have a simple and efficient method. The 
stencil is narrower than for the B-spline formulas that were given previously. 
With the formula above we also get a symmetric stencil, which is more natural if 
the other parts of the computation, such as difference approximations of PDEs 
are done by symmetric formulas. Furthermore, symmetry makes periodic BCs 
somewhat easier to implement. Note that the absence of symmetry for either 
the scaling function or the wavelet can lead to phase distortion. This can be 
shown to be important in signal processing applications. 


Multi-Dimensional Wavelets. The computation of multi-dimensional wavelets 

is quite expensive, especially in 3-D. A simple minded efficient way is to eval- 
uate the wavelet coefficients dimension-by-dimension. This means that we get 
two set of wavelet coefficients w^ ^y) and w y m k {x), where now (j, k) is the 
position and m is the scale. The precise definition is 

w mj(y) = f f(x,y)lpm,j{x)dx 

w m,k( x ) = I f ( x > y)^m,k{y) d y- 

Thus, the dimension-by-dimension approach involved only terms evaluated 
as finite differences in the ^-direction and terms which are evaluated in the y- 
direction. We then use the ( y ) coefficients for the x-direction computation, 

and the y-coefficients for the y-direction computation. 

Shock/Shear Wavelet Sensor. For the numerical experiments presented in 
the next section the wavelet sensor is obtained by computing a vector of the 
approximated Lipschitz exponent of a chosen vector function to be sensed 
with a suitable multiresolution non-orthogonal wavelet basis function. Here, 
“vectors or variables to be sensed” means the represented vectors or variables 
that are suitable for the extraction of the desired flow physics. The variables 
to be sensed can be the density, the combination of density and pressure, the 
characteristic variables, the jumps in the characteristic variables 5*. +1 y 2 , or the 
entropy variable vector W ([20, 83]). 
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For example, if the characteristic variables are the chosen vector to be sensed 
by the wavelet approach, the sensor can be defined as 

+ 1/2 ~ r ( a j+ 1 / 2)5 (5.26) 

where a*- +1 / 2 estimated Lipschitz exponent of the Ith characteristic com- 
ponent with l — 1,2, 3,4, the four characteristic waves. r(a) is a sensing 
function which decreases from r(0) = 1 to r(l) = 0 (for the aforementioned 
type of wavelets). Note that the Zth component of the estimated Lipschitz expo- 
nent a l j+i / 2 is not to be confused with the jump in the Zth characteristic variables 

ctj+i / 2 in Section 5.1. We use <Sj +1 / 2 as the sensor to distinguish it from the 
ACM sensor u/. +1 / 2 * n Section 5.1. Noted that the k index is omitted (for the 
2-D case) for simplicity. 

If we instead base the exponent estimate on point centered quantities, we 
will use the sensor function 


S l j+ 1/2 = max(r(o5-),r(a5 +1 )). (5.27) 

If the exponent estimate is based on other quantities than the characteristic 
variables, (e.g., density and pressure), we use the switch 

1/2 — max<Sj + 1 / 2 ? (5.28) 


where the maximum is taken over all components of the waves used in the 
estimate. In this case, the switch is the same for all characteristic fields. 

The function r(a) should be such that r(0) = 1, and r(l) = 0. Three 
options considered are 


T 



1 

0 


a < ao 
a > ao 


r(a) = ^ ^ arctan K(oto — ot) 


(5.29) 


r (a) = max{0, min[l, (a — 1)/ (ao — 1)]}. 


Here, ao is a cut off exponent to be chosen. For the arctan function the values 
0 and 1 are not attained, but we take the constant K large enough so that the 
function is close to zero for a > 1, and close to one for a < 0. We have 
tried values for K in the interval [200,500]. Alternatively, one can integrate 
the actual a value into the sensor function instead of using the same amount of 
numerical dissipation at the cut off exponent. 

After some experimentation we have found that switching on the dissipation 
at the grid points where a < 0.5 works well, i.e., 


r(a) = 


1 a < 0.5 

0 a > 0.5 


(5.30) 
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In fact the method does not seem to be very sensitive to the exact value of cut 
off ao, (for 0.4 < a 0 < 0.6) for all the test cases considered. Furthermore, the 
same cut off value, 0.5, works well for all problems we have tried in Section 6 
(except for the vortex convection case, where ao = 0.0 is used in conjunction 
with entropy splitting [83]). Experiments with smoothed step functions do not 
give very different results. 

5.3 Test Example of the Shock/Shear Wavelet Detectors 

This section shows the performance of the wavelet sensor using the dimension- 
by-dimension approach for a 2-D complex flow structure. It is important to note 
that the illustration involves only the feature extraction capability of the wavelet 
sensor on a given grid function. No dynamic behavior was involved (i.e., the 
numerical scheme is not part of the analysis). Figure 5.1 shows the computed 
density and pressure contours from a precomputed numerical simulation at 
t = 120 with At — 0.12 to be used as the two-dimensional discrete functions 
to be analyzed by the wavelet algorithm. The discrete functions represent a 
numerical data of a shock from the upper left comer, impinging on a horizontal 
shear layer in the middle of the domain. The shock is reflected from the lower 
wall boundary. For more details about the problem, see Yee et al. [82, 83]. 

Density 



Fig. 5.1. 2-D Testing discrete function, (density and pressure contours at 
t = 120). 
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Figure 5.2 shows contours of the estimated Lipschitz exponent a for the function 
in Fig. 5.1. The value a was computed here from three levels (mo = 3) of the 
wavelet algorithm, using the wavelet coefficient 

w m,j,k = y/(Kij,k ) 2 + t) 2 ’ < 5 - 31 ) 

where the one dimensional coefficients were computed by the multiresolution 
operators (5.24) in each coordinate direction. The coefficients were computed 
for the pressure. The top figure in Fig. 5.2 shows a contours on levels from 
0.5 to 0.9. The lower figure shows the corresponding sensor, a function which 
is one for a < 0.5 and zero otherwise. The wavelet sensor clearly captures 
the shock and the shear layer. The low a at the upper boundary to the right is 
probably due to mildly unstable BCs at the upper boundary. 


Alpha contours [0.5 0.9] 



Sensor, contour at 0.5 



Fig. 5.2. Top: a contours 0.5 < a < 0.9; Bottom: sensor contour at 
a = 0.5. by the RH -wavelet. 


5.4 Blending of Different Filters 

The nonlinear filters for the ACM or wavelet shock/shear capturing nonlinear 
filter might not be sufficient for (a) time-marching to steady state and (b) spuri- 
ous high frequency oscillations due to insufficient grid resolution and nonlinear 
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instability away from discontinuities, especially for turbulent and large-eddy 
simulations. This Section discusses the blending of other filters with these 
shock/shear filters. 

In classical CFD codes, a second order accurate base method is used together 
with two constant strength linear numerical dissipation terms. One linear fourth- 
order dissipation is used everywhere except near shocks/shears/steep-gradients 
to remove nonlinear instabilities. It does not affect the second order accu- 
racy of the base scheme. The second dissipation term is a second-order linear 
dissipation, which affects the order of accuracy, but is only switched on near 
discontinuities, and/or steep unresolved gradients using a gradient sensor. The 
sensor used cannot distinguish the different flow features distinctly and is not 
accurate enough for turbulent statistics and long-time acoustic computations, 
unless extreme grid refinement is employed. 

In analogy with the aforementioned classical methods, a more advanced 
numerical dissipation model with improved flow feature extraction sensors for 
high order central schemes is proposed. Here, we consider a dissipation model 
with two parts. One part is a nonlinear filter ([82]) and the second part is a 
high order linear numerical dissipation term modified at boundaries to become 
a semi-bounded operator, see [67, 65]. The wavelet dissipation control sensor 
developed in [63] is used as the flow feature detector. 

Time-Marching to Steady State. For time-marching to steady state one 
usually needs to add fourth-order linear dissipation to a second-order spatial 
differencing scheme (Beam and Warming (1976)). For the present schemes us- 
ing characteristic filters, in addition to the nonlinear shock/shear filter, one might 
need to add a sixth-order linear dissipation to a fourth-order spatial base scheme 
and an eighth-order linear dissipation to a sixth-order spatial base scheme in 
regions away from shocks for stability and convergence. Let L d be such an 
additional filter operator. The two ways of incorporating the L & operator are 
options I and II discussed in Section 4. 

To minimize the amount of dissipation due to L d in the vicinity of shock 
waves, there should be a switching mechanism K& to turn off the L d operator in 
the vicinity of shock waves. The Lj operator can be applied to the conservative, 
primitive or characteristic variables. The simplest form is to apply to the 
conservative variables. Alternatively, since all of the work in computing the 
average states and the characteristic variables is done for the shock-capturing 
filter operator, one can apply the L d operator to the characteristic variables. 
In this case, the switching mechanism k ^ can be a vector so that it is more in 
tune with the nonlinear shock detector using the approximate Riemann solver. 
For example, one can set k = 0 for the linearly degenerate fields and blend a 
small amount of to remove spurious noise generated by the lack of nonlinear 


Adaptive Low-Dissipative Schemes 


31 


filters. This blending of the nonlinear shock/shear filter with the L ^ operator 
can be applied to time-accurate computations as well. 

Suppression of Spurious High Frequency Oscillations. The nonlinear 
shock/shear filters might not be able to remove spurious high frequency oscilla- 
tions effectively unless sufficient fine grid points are used. For the suppression 
of unphysical high frequency oscillations due to insufficient grid resolution 
and nonlinear instability away from discontinuities, higher-order spectral-like 
filters (Vichnevetsky (1974), Lele (1994), Alpert (1981), Visbal and Gaitonde 
(1998), Gaitonde and Visbal (1999)) might be needed at the locations where 
the value of the shock/shear sensor is very small or zero. If spectral-like filters 
are needed, a proper blending of nonlinear shock/shear filters with spectral-like 
filters should be applied. In this case, we can use the same procedures as the 
time-marching to the steady state except the La operator should be replaced 
with the spectral-like filters (for compact central schemes). 


An Adaptive Numerical Dissipation Model for Shock-Turbulence Inter- 
actions. Assume a sufficient grid is used for the problem and scheme in 
question, and that the scale of turbulent fluctuations is larger than the spuri- 
ous high-frequency oscillations. Below we present a filter model under these 
assumptions. If the scale of the turbulent fluctuation is in similar scale as the 
high-frequency oscillations, a different wavelet with a turbulent fluctuation sen- 
sor should be added. For example, for a sixth-order central spatial base scheme, 
we define the 1-D filter numerical flux of the numerical dissipation operator as 

Li d 

H j- 1/2- 

Hf_ 1/2 = S,_ 1/2 F;_ 1/2 + dj[ 1 - Sj_ l/2 ](h 6 D-(D + D-) 3 Uj, (5.32) 

here Sj_i /2 is a switch computed as described in Section 5.2.4, and Fj _\/2 > s 
the flux function corresponding to the dissipative portion of a shock-capturing 
scheme(e.g., second order accurate TVD scheme) [82]. The first part of the filter 
stabilizes the scheme at shock/shear locations. The second part is an eighth- 
order linear filter which improves nonlinear stability away from shock/shear 
locations. Analogous eighth-order filters can be used if a sixth-order compact 
spatial base scheme is used [18, 76]. We switch on the high order part of the 
filter when we switch off the nonlinear filter. The physical quantity (e.g., local 
Mach number) can be used to determine the d,j parameter of this high order 
dissipation term. 

To further increase stability properties, it is possible to use the sensor to switch 
on and off the entropy splitting and adjust the value of the entropy splitting 
parameter according to flow type and region. For the 1-D shock/turbulence 
interactions to be presented in the next section, however, we believe a constant 
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/? = 1 away from the shock waves is sufficient. After the completion of a full 
time step computation using the sixth-order base scheme (denoting the solution 
by Uj), we filter this solution by 


v T 1 = °i + T [fl i+>/2 “ (533) 


Here the filter numerical fluxes Hj ±1 ^ 2 are evaluated at U. 

5.5 Spurious Numerical Solutions and Instability Due to 
Under-Resolved Grids 

There has been much discussion on verification and validation processes for 
establishing the credibility of CFD simulations [68, 7, 75, 21, 52]. Since the 
early 1990s, many of the aeronautical and mechanical engineering related ref- 
erence journals mandated that any accepted articles in numerical simulations 
(without known solutions to compare with) need to perform a minimum of one 
level of grid refinement and time step reduction. On the other hand, it has 
become common to regard high order schemes as more accurate, reliable and 
requiring less grid points. The danger comes when one tries to perform com- 
putations with the coarsest grid possible while still hoping to maintain numeri- 
cal results sufficiently accurate for complex flows and, especially, data-limited 
problems. On one hand, high order methods when applied to highly coupled 
multidimensional complex nonlinear problems might have different stability, 
convergence and reliability behavior than their well studied low order counter- 
parts, especially for nonlinear schemes such as TVD, MUSCL with limiters, 
ENO, WENO, and spectral elements and discrete Galerkin. See for example 
references [23, 74, 6, 49, 78, 81, 80, 79]. On the other hand, high order meth- 
ods involve higher operation counts per grid and systematic grid convergence 
studies can be time consuming and prohibitively expensive. At the same time it 
is difficult to fully understand or categorize the different nonlinear behavior of 
finite discretizations, especially at the limits of under-resolution when different 
types of numerical (spurious) bifurcation phenomena might occur, depending 
on the combination of grid spacings, time steps, initial conditions (ICs) and 
numerical treatments of BCs as well as the nonlinear stability of the scheme in 
question. 
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Due to the difficulty in analysis, the effect of under-resolved grids and the 
nonlinear behavior of available spatial discretizations are scarcely discussed 
in the literature. Here, an under-resolved numerical simulation, according to 
Brown & Minion, is one where the grid spacing being used is too coarse to re- 
solve the smallest physically relevant scales of the chosen continuum governing 
equations that are of interest to the numerical modeler. Before the nineties, it 
was common in DNS to avoid the use of numerical dissipations. It was standard 
practice to refine the grid not just to resolve the multiscale physics but also to 
overcome nonlinear instability instead of employing numerical dissipation or 
filters. In certain cases, the grid is finer than what is needed to resolve the 
smallest scale. This section illustrates the situation where the necessity of finer 
grid can be overcome by the use of an appropriate filter and still be able to 
obtain an accurate and stable solution with a much coarser grid. 

Brown and Minion [6, 49] studied the effects of under-resolved grids by 
considering the shear layer roll-up problem that arises when the Navier-Stok.es 
equations are solved in the unit square with doubly-periodic BCs with ICs given 
by 

_ f tanh(p(y - 0.25)) for y < 0.5 (5.34) 

U — \ tanh(p(0.75 — y )) for y > 0.5 

v — 0.05 sin(27rx). (5.35) 

In [6, 49], the behavior of several difference methods was considered. These 
difference methods include a Godunov projection method, a primitive variable 
ENO method, an upwind vorticity stream-function method, centered difference 
methods of both a pressure-Poisson and vorticity stream-function formulation, 
and a pseudospectral method. They demonstrated that all these methods pro- 
duce spurious, non-physical vortices. While these extra vortices might appear 
to be physically reasonable, they disappear when the mesh is refined. 

Figure 5.3 shows filter-based spectral element results for the problem (5.34) 
as computed by Fischer and Mullen [17]. The spectral element method is char- 
acterized by the discretization pair (E, N), where E is the number of quadri- 
lateral elements and N is the order of the tensor-product polynomial expansion 
within each element. This filter presented in [17] is designed to stabilize the 
PN2 spectral element method at high Reynolds numbers. The PN 2 method, 
introduced by Maday and Patera [40] , is a consistent approximation to the Stokes 
problem which employs continuous velocity expansions of order N and discon- 
tinuous pressure expansions of order N - 2. The discretizations in Fig. 5.3a - 
5.3e consist of a 16 array of elements, while Fig. 5. 3f consists of a 32 x 32 array. 
Here, a denotes the spectral filter coefficient (not to confused with the Lipschitz 
exponent or the jump in the characteristics in Sections 5.1 and 5.2), with a — 0 
corresponding to no filtering. The time step size is At = 0.002 in all cases, 
corresponding to CFL numbers in the range of 1 to 5. Without filtering, Fischer 
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and Mullen were not able to simulate this problem at any reasonable resolution. 
Figure 5.3a illustrates the result prior to blow up for the unfiltered case with 
(E,N) = (16 2 , 16), which has a resolution corresponding to an n x ngrid with 
n = 256. Unfiltered results for N = 8 (n = 128) and TV = 32 (n = 512) are 
similar. Filtering with a — 0.3 yields dramatic improvement for n = 256 (Fig. 
5.3b) and n — 128 (Fig. 5.3d). Though the so-called full projection with filter 
strength a = 1 is stable, the partial filtering of (a < 1) gives smoother results 
and is preferable. The cases in 5.3e and 5.3f correspond to the difficult “thin” 
shear layer case of [6] and show the benefits of high-order discretizations. Both 
cases correspond to a resolution of n 2 = 256 2 . In Fig. 5.3e, this is attained 
with (Ej N ) = (16 2 , 16), while in Fig 5.3f, (£*, N) — (32 2 , 8). Although both 
results are stable (due to the filter), Fig 5.3f reveals the presence of spurious 
vortices that are absent in the higher-order case. 

6. Numerical Examples 

This section illustrates the power of entropy splitting, the difference in per- 
formance of linear and nonlinear (with sensor controls) filters and the combi- 
nation of both types of filters with adaptive sensor controls. We use the same 
notation as in [82, 83, 64]. The artificial compression method (ACM) and 
wavelet filter schemes using a second-order nonlinear filter with sixth-order 
spatial central interior scheme for both the inviscid and viscous flux derivatives 
are denoted by ACM66 and WAV66. See [82, 83, 64] for the forms of these 
filter schemes. The same scheme without filters is denoted by CEN66. The 
scheme using the fifth-order WENO for the inviscid flux derivatives and sixth- 
order central for viscous flux derivatives is denoted by WEN05. Computations 
using the standard fourth-order Runge-Kutta temporal discretization are indi- 
cated by appending the letters “RK4” as in CEN66-RK4. ACM66 and WAV66 
use the Roe’s average state and the van Leer limiter for the nonlinear numer- 
ical dissipation portion of the filter. The wavelet decomposition is applied in 
density and pressure, and the maximum wavelet coefficient of the two com- 
ponents is used. The nonlinear numerical dissipation is switched on wherever 
the wavelet analysis gives a Lipschitz exponent [63] less than 0.5. Increasing 
this number will reduce oscillations, at the price of reduced accuracy (see [63] 
for other possibilities). Computations using entropy splitting are indicated by 
appending the letters “ENT” as in WAV66-RK4-ENT. Computations using an 
eighth-order linear dissipation filter are indicated by appending the letters “D8” 
as in WAV66-RK4-D8. In order not to introduce additional notation, inviscid 
flow simulations are designated by the same notation, with the viscous terms 
not activated. 
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Isentropic Vortex Evolution 

(Horizontally Convecting Vortex, vortex strength 0=5) 


Freestream : 

(tioo^oo) = (1,0); Poo Poo 1 
JC: Perturbations are added to the freestream (not in entropy) 

( Su,Sv ) - A- e — (-(y - y 0 ),(x - * 0 )) 

ST = 1 )^ e l-r» 

87 5T 2 

1 2 = (x-x/ + (y-yo? 

Computational Domain & Grid Size : 


0 < * <10 & -5 < y < 5 
80 x 79 Uniform grid 
Periodic BC in x & y 


Initial Vortex, Density Contours 



( 2 ®o> y*o) 
center of vortex 


Figure 6. 1. Vortex convection problem description. 


6.1 A 2-D Vortex Convection Model [82, 83, 63, 65] 

The onset of nonlinear instability of long-time numerical integration, the 
benefit of the entropy splitting and the difference in performance of linear and 
nonlinear numerical dissipations in improving nonlinear stability for a horizon- 
tally convecting vortex (see Fig. 6.1) can be found in in [82, 83, 63, 65]. We 
summarize the results here. 

To show the onset of nonlinear instability, the 2-D perfect gas compressible 
Euler equations are approximated by CEN66-RK4 with periodic BCs imposed 
using a 80 x 79 grid with the time step At = 0.01. Since this is a pure 
convection problem, the vortex should convect without any distortion if the 
numerical scheme is highly accurate and non-dissipative. Although CEN66- 
RK4 is linearly stable, the test problem is nonlinear and instability eventually 
sets in. Almost perfect vortex preservation is observed for up to 5 periods of 
integrations (5 times after the vortex has convected back to the same position - 
time = 50). Beyond 5 periods the solution becomes oscillatory, and blows up 
before the completion of 6 periods. The blow up is associated with an increase 
in entropy [65]. If we instead use the entropy-split form of the approximation 
(CEN 66-RK4-ENT) with a split parameter (3=1, almost perfect vortex preser- 
vation for up to 40 periods can be obtained. The computation remains stable for 
up to 67 periods before it breaks down. The time history of the L 2 entropy norm 
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Entropy Norm vs Time, CEN66-RK4-ENT [3=1 



Figure 6.2. Entropy norm history of CEN66-ENT: entropy split parameter (3 — 1 and 80 x 79 
grid. 


and density contours of the IC and the computed solution after 5, 10, 30, 50 and 
67 periods using CEN66-RK4-ENT is shown in Figs. 6.2 and 6.3. The norm 
is decreasing, although the instabilities break down the solution after 67 peri- 
ods. Using the second-order nonlinear filter without splitting (ACM66-RK4 
or WAV66-RK4), the solution remains stable beyond 67 periods. However, 
the numerical solution gradually starts to diffuse after 20 periods. If we use 
the nonlinear filter in conjunction with entropy splitting (ACM66-RK4-ENT or 
WAV66-RK4-ENT), almost perfect vortex preservation can be obtained for up 
to 120 periods using a split parameter (3=1 [83]. 

The density contours of the solution after 5, 10, 200 and 300 periods for 
the un-split (/3 = oo) computation using the eighth-order linear dissipation 
(CEN66-RK4-D8) are shown in Fig. 6.4. The linear dissipation {-dh 7 (D+D^Uj) 
with grid spacing h was added to the sixth-order base scheme to discretize the 
convection terms. The parameter d is a given constant (d = 0.002) and is 
scaled with the spectral radius of the Jacobian of the flux function, and D + 
and -D- are the forward and backward difference operators, respectively. This 
numerical dissipation is applied as part of the scheme and not as a post pro- 
cessing filter. The computation can be run for 300 periods without breakdown. 
However, serious degradation of accuracy occurs after 250 periods. For this 
particular problem, the CEN66-RK4-D8 out performed the ACM66-RK4-ENT 
and WAV66-RK4-ENT using (3=1 . Perhaps using a higher than third-order 
nonlinear filter might improve the performance of the ACM66-RK4-ENT and 
WAV66-RK4-ENT at the expense of an increase in CPU. 
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Figure 6.5. Mean velocity profiles 


6.2 DNS of 3-D Compressible Turbulent Channel Flow 

[61] 

This numerical example illustrates the power of entropy splitting for DNS 
computations. To obtain accurate turbulent statistics, very long-time integration 
and highly accurate methods are required for this DNS computation. The com- 
putation employed the SBP-satisfying boundary difference operator [9] with 
the fourth-order central interior scheme applied to the split form of the inviscid 
flux derivatives CEN44-RK4-ENT with /3 = 4, and aLaplacian viscous formu- 
lation. The fluid mechanics of this 3-D wall bounded isothermal compressible 
turbulent channel flow has been studied in some detail by Coleman et al. [11]. 
They showed that the only compressibility effect at moderate Mach numbers 
comes from the variation of fluid properties with temperature. They used a uni- 
form body force term to drive the flow, but recommended the constant pressure 
gradient approach which was adapted by Sandham et al. [61]. 

Grid refinement study. A simplified case is taken in which the fluid proper- 
ties (viscosity and thermal conductivity) are held constant and the computational 
box size is kept small. The latter is justified as a method of reducing cost as the 
gross turbulence statistics are relatively insensitive to the computation box size, 
so long as the domains are still significantly larger than the minimal domains on 
which turbulence can be sustained. A Mach number of 0. 1 is chosen, based on 
friction velocity and sound speed corresponding to the fixed wall temperature. 
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Figure 6.6. Effect of grid refinement on normal stresses. The top curves relate to the left scale, 

the middle to the right scale and the lowest to the furthest right scale. 

Channel half width h, friction velocity u T , wall temperature and bulk (inte- 
grated) density are the normalizing quantities for non-dimensionalization with 
a Reynolds number of 180. Together with the constant property assumption, 
this choice of Mach number means that results can reasonably be compared 
to results from previous incompressible flow calculations. The computations 
were carried out at a fixed CFL=2.0. They were started with artificial ICs and 
first run to time t = 50, by which time dependence on the ICs is lost. Statistics 
were accumulated over the time interval t = 200 to t — 300. 

Three grids 12 x 41 x 12, 24 x 81 x 24 and 36 x 121 x 36 were considered. 
The largest number in each case corresponds to the direction normal to the wall 
(■ y ). The computational box has non-dimensional length 3 in the (streamwise) in- 
direction, 1.5 in the (spanwise) z-direction and 2 in the y-direction. The x- and 
z-directions have periodic BCs with uniform grid spacing. In the (/-direction, 
the grid is stretched according to 

y _ tanh(c t? 7?) 

h tanh 

with r] uniformly distributed on [-1,1], c rj = 1.7. The ratio of grid points in 
each direction was chosen so that all directions have roughly the same degree 
of resolution of the relevant turbulence microscales in each direction. Figure 
6.5 shows the mean flow velocity, Fig. 6.6 the root mean normal stresses and 
Fig. 6.7 the stress profiles across the channel. Angle brackets () denote averages 
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Figure 6.7. Effect of grid refinement on turbulent shear stress (curve falling to zero at the 
walls). The total stress (straight lines, non-zero at walls) is also shown. 


over the homogeneous spatial directions and time while in the usual notation 
double primes denote deviations from mass-weighted (Favre) averages. The 
convergence is not uniform across the channel but the change from medium to 
fine grid is smaller than the change from coarse to medium grid. A comparison 
of the rms quantities with an incompressible flow simulation on the same size 
computational box (Z. Hu, private communication) is shown in Fig. 6.8. Here 
we compare the 36 x 121 x 36 fourth-order compressible simulation with a 
32 x 81 x 32 fully spectral incompressible simulation using the method described 
in [59]. As expected, good agreement is found, as expected for this Mach 
number (0.1 based on friction velocity or 1.8 based on centerline velocity). 

The convergence of various global measures can be found in Table 1 of [61] 
for the three grids. For the pressure gradient and Reynolds number specified, the 
velocity gradient at the wall should be 180, the difference away from this being 
an error of the simulation. Here Re r is the Reynolds number based on u T , the 
mean density at the wall (p w ) and the mean viscosity (fi w ) at the wall. For the 
finest grid the resolution in wall units (a common check on resolution in DNS) 
is A+ = 15 and A+ = 7.5 and approximately 10 points are in the sublayer 
y+ < 10. The simulations demonstrate a robustness down to very coarse 
resolutions, comparable with the best incompressible turbulent flow solvers 
incorporating de-aliasing and skew-symmetric formulation of the convective 
terms. Without the use of the entropy splitting of the inviscid flux derivatives 
and without the use of the Laplacian viscous formulation, the CEN44-RK4 
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Figure 6.8. Comparison of weighted turbulence quantities on a 36 x 121 x 36 grid with an 
incompressible flow simulation on a 32 x 32 x 81 grid. 



Figure 6.9. Mean velocity profile, comparing current simulation (dashed line) with Coleman 
et < 2 /.[ll](solid line). 


(un-split) solutions for the same CFL number, diverge for all three grids before 
meaningful turbulence statistics can be obtained. 
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Figure 6.10. Root mean normal turbulent stresses, comparing current simulation (dashed line) 
with Coleman et al. [ 1 1 ](solid line). The top curves relate to the left scale, the middle to the right 
scale and the lowest to the furthest right scale. 


Comparison with Coleman et al.. Coleman et al [11] carried out compa- 
rable simulations in their study of the effects of Mach number on turbulence 
statistics. This section shows the simulation of their case Re T = 190 and 
M t = 0.095 with a uniform body force term together with variable fluid prop- 
erties (power-law temperature dependence of the viscosity with exponent 0.7 
and fixed Prandtl number Pr — 0.7). With the variable viscosity there is a need 
to use a larger computational box size than was used in the preceding section, 
since turbulence structures become larger as the viscosity is reduced (the wall is 
cold relative to the bulk flow). We chose to use a box of size 6x2x3, i.e., twice 
as large in x and z as in the preceding section. This size is still somewhat lower 
than that of Coleman et al . , who used a box of size 47r x 2 x 47 t/3. A computa- 
tional grid of 60 x 141 x 60 was used, giving Ax + = 19 and A z + — 9.5. These 
are comparable to those used by Coleman et al (16.6 and 10.0 respectively). 
There were 12 points in the sublayer (y+ < 10). 

For this simulation a parallel implementation was used, which illustrated the 
excellent parallel scaling of the method on a Cray T3E-1200E computer (90% 
efficiency for a 240 3 benchmark on 256 processors and continued good scaling 
up to 768 processors, as reported in Ashworth et al [3]). The simulation 
presented here used 32 processing elements. 

Table 1 shows a summary of the output from the simulation. Data from 
Coleman et al have been re-normalized for comparison with the current simu- 
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Figure 6.1 1. Turbulent and total shear stresses, comparing current simulation (dashed line) 
with Coleman et al[ 1 l](solid line). 


Table 1. Comparison of centerline velocity U c > bulk velocity Ub, wall shear stress (p^pj)w and 
shape factor H with Coleman et al. [1 1] 


Simulation 

Uc 

u b 

) w 

H 

Current 

18.9 

16.3 

190.3 

1.66 

Coleman et al. 

18.5 

15.9 

189.5 

1.65 


lation results. Figures 6.10 and 6.11 show the shear stress and rms turbulence 
fluctuations, while Fig. 6.9 shows the mean velocity profile. Overall a good 
agreement is obtained illustrating the good performance of the method for a 
resolution comparable to that of a spectral method. Good turbulence kinetic 
energy budgets have also been obtained [38]. 

We note that for this well-studied problem with accurate turbulent flow 
databases for comparison, we can safely conclude that entropy splitting in con- 
junction with the Laplacian formulation calculations was able to obtain stable 
and fairly accurate solutions using coarse to moderate grid sizes without added 
numerical dissipation or filters. Unlike the spectral method, this high order 
method can be efficiently extended to general geometries [77]. For the same 
3-D problem, the finite difference formulation of the WEN05 is more than six 
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times as expensive, yet more diffusive than the present scheme using the same 
temporal discretization. 

The numerical methods are currently being applied to several practical prob- 
lems. Alam & Sandham [1] have studied shock-free transition to turbulence 
near the leading edge of an aerofoil, while Lawal & Sandham [35] have used 
the above method in conjunction with the high order nonlinear filter scheme 
from Yee et cil. [82] to study transitional shock boundary layer interaction in 
flow over a bump. These practical applications have been run without the need 
for changes to the numerical method and hence are leading to some confidence 
that the developments presented here are generally applicable for DNS of com- 
pressible turbulent flow . 

For the performance of ACM66-RK4, WAV66-RK4 and WEN05-RK4 on a 
spatially or a time-developing mixing layer problem containing shock waves, 
see [83, 63]. 

6.3 Computational Aeroacoustics Applications (CAA) 

[50] 

This numerical example illustrates the applicability of the entropy splitting 
to CAA for low Mach number flows. The numerical prediction of vortex sound 
has been an important goal in CAA since the noise in turbulent flow is generated 
by vortices. To verify our numerical approach for CAA, the Kirchhoff vortex 
is chosen for the numerical test. The Kirchhoff vortex is an elliptical patch of 
constant vorticity rotating with constant angular frequency in irrotational flow. 
The acoustic pressure generated by the Kirchhoff vortex is governed by the 
2D Helmholtz equation, which can be solved analytically for almost circular 
Kirchhoff vortices using separation of variables. See [50] and references cited 
therein for details. The difficulty with this test case is the large gradient of the 
acoustic pressure adjacent to the Kirchhoff vortex. 

The perturbation form of the entropy split 2D Euler equations in conjunction 
with a fourth-order linear filter operator CEN66-RK4-ENT-D4 was applied to 
Kirchhoff vortex sound at low Mach number using a high order metric evalua- 
tion of the coordinate transformation. SBP operators using the SAT method of 
implementing the time-dependent physical BC was employed. Due to the large 
disparity of acoustic and stagnation quantities in low Mach number aeroacous- 
tics, the split Euler equations are formulated in perturbation form to minimize 
numerical cancellation errors. 

A very accurate numerical solution with a relatively coarse grid was ob- 
tained using CEN66-RK4-ENT-D4 compared with the un-split (CEN66-RK4- 
D4), and un-filter cases CEN66-RK4, CEN22-RK4 and CEN44-RK4. The 
extra CPU due to the use of the split form of the inviscid flux derivatives is 
more than compensated by the improved accuracy and stability of the numer- 
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ical simulation, especially near regions of large gradients. For this weakly 
nonlinear test case that does not require long time integration, the amount of 
filter needed, although very small, is still important. Higher order filter opera- 
tors and the formulation of the numerical BCs and filter operators in terms of 
the entropy variables to satisfy a discrete energy estimate in a nonlinear sense 
will be considered in the future. 

6.4 Multiscale Complex Unsteady Viscous Compressible 
Flows [64, 66] 

Extensive grid convergence studies using WAV 66-RK4 and ACM66-RK4 for 
two complex highly unsteady viscous compressible flows are given in [64, 66]. 
The first flow is a 2-D complex viscous shock/shear/boundary-layer interaction. 
This is the same problem and flow conditions studied in Daru & Tenaud [12]. 
The second flow is a supersonic viscous reacting flow concerning fuel breakup. 
More accurate solutions were obtained with WAV66-RK4 and ACM66-RK4 
than with WEN05-RK4, which is nearly three times as expensive. To illustrate 
the performance of these nonlinear filter schemes, the first model is considered. 
The ideal gas compressible full Navier-Stokes equations with no slip BCs at 
the adiabatic walls are used. The fluid is at rest in a 2-D box 0 < x,y < 1. 
A membrane with an initial shock Mach number of 2.37 located at x = 1/2 
separates two different states of the gas. The dimensionless initial states are 

pi = 120, PI , = 120/7; p R = 1.2, pr = I.2/7, (6-1) 

where pl,Pl are the density and pressure respectively, to the left of a: = 1/2, 
and pr, pr are the same quantities to the right of x = 1/2. 7 = 1.4 and the 
Prandtl number is 0.73. The viscosity is assumed to be constant and independent 
of temperature, so Sutherland’s law is not used. The velocities and the normal 
derivative of the temperature at the boundaries are set equal to zero. This is 
done by leaving the value of the density obtained by the one sided difference 
scheme at the boundary unchanged, and updating the energy at the boundary 
to make the temperature derivative equal to zero. 

At time zero the membrane is removed and wave interaction occurs. An 
expansion wave and a shock are formed initially. Then, a boundary layer is 
formed on the lower boundary behind the right going waves. After reflection, the 
left going shock wave interacts with the newly formed boundary layer, causing 
a number of vortices and lambda shocks near the boundary layer. Other kinds of 
layers remain after the shock reflection near the right wall. The complexity of 
this highly unsteady shock/shear/boundary-layer interactions increases as the 
Reynolds number increases. 

As an illustration, we show here the difficult case of Reynolds number 
Re = 1000. The computations stop at the dimensionless time 1 when the 
reflected shock wave has almost reached the middle of the domain, x = 1/2. 
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The numerical results discussed here are at time 1 with a uniform Cartesian grid 
spacing as described by Daru and Tenaud. Due to symmetry, only the lower 
half of the domain is used in the computations; symmetry BCs are enforced 
at the boundary y — 1/2. Figure 6.12 shows the comparison of a second- 
order MUSCL using a second-order Runge-Kutta method (MUSCL-RK2) with 
WAV66-RK4, ACM66-RK4 and WEN05-RK4 using a 1000 x 500 grid. Com- 
paring with the converged solution of WAV66-RK4 and ACM66-RK4 using 
3000 x 1500 (see bottom of figure) and 4000 x 2000 grids (see [64]), one 
can conclude that WAV66-RK4 exhibits the most accurate result among the 
1000 x 500 grid computations. We note that, for this Reynolds number, the 
unsteady problem is extremely stiff, requiring very small time steps and very 
long-time integrations before reaching the dimensionless time of 1. 

6.5 1-D Shock-Turbulence Interactions Using the 

Adaptive Numerical Dissipation Model 

The dissipative model (5.32) is used to solve a simple, yet difficult, 1-D 
compressible inviscid shock-turbulence interaction problem with initial data 
consisting of a shock propagating into an oscillatory density. The initial data 
are given by 


{p L , u L , p L ) = (3.857143, 2.629369, 10.33333) (6.2) 

to the left of a shock located at x — —4, and 

( PR , ur, Pr) = (1 + 0.2sin(5:r), 0, 1) (6.3) 

to the right of the shock where u is the velocity. Fig. 6.13 show the compar- 
ison between a second-order MUSCL-RK2 with a sixth-order central scheme 
and the aforementioned numerical dissipation model using RK4 as the time 
discretization (WAV66-RK4-D8). The parameter d = 0.002 is scaled with the 
spectral radius of the Jacobian of the flux function. Note that the eighth-order 
dissipation is a filter, and is different from the CEN66-D8 used in Section 6.3 
where the dissipation is part of the scheme. The solution using a second-order 
uniformly non-oscillatory (UNO) scheme on a 4000 uniform grid is used as 
the reference solution (solid lines on the first three sub-figures). The bottom 
of the right figures show the density and Lipschitz exponent distribution for 
the WAV66-RK4-D8 using 400 grid points. Comparing our result with the 
most accurate computation found in the literature for this problem, the current 
approach is highly efficient and accurate using only 800 grid points without 
grid adaptation or a very high order shock-capturing scheme. For the present 
computation, the WAV66-RK4-D8 consumed only slightly more CPU than the 
second-order scheme MUSCL-RK2. With the eighth-order dissipation filter 



Density at time 1 



(a) MUSCL-RK2, 1000 x 500 grid 


Density at time 1 
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(c) WAV66-RK4, 1000 x 500 grid 


WAV66 3000x1500 



(e) WAV66-RK4, 3000 x 1500 grid 



(f) ACM66-RK4, 3000 x 1500 grid 


Figure 6.12. Comparison: MUSCL-RK2, WAV66-RK4, WAV66-RK4 and WEN05-RK4 for 
Re — 1000. Density contours using 1000 x 500 and 3000 x 1500 grids. 
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(c) WAV66-RK4-D8, 800 grid 
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(d) Density & Lipschitz Exponent Distributions, 
WAV66-RK4-D8, 400 grid 


Figure 6.13. Comparison and Lipschitz exponent distribution of WAV66-RK4-D8, second- 
order B -spline wavelet. 

turned off (i.e., only the nonlinear filter remains - WAV66-RK4), the computa- 
tion is not very stable unless a finer grid and smaller time step is used. Turning 
on the entropy splitting away from the shocks helps to reduce the amount of the 
eighth-order dissipation coefficient [65]. 

7. Concluding Remarks 

An integrated approach for the control of the numerical-dissipation/filter 
in high order schemes for numerical simulations of multiscale complex flow 
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problems is presented. The approach is an attempt to further improve nonlinear 
stability, accuracy and efficiency of long-time numerical integration of com- 
plex shock/turbulence/acoustics interactions and numerical combustion. The 
required type and amount of numerical-dissipation/filter for these flow prob- 
lems are not only physical problem dependent, but also vary from one flow 
region to another. Among other design criteria, the key idea consists of au- 
tomatic detection of different flow features as distinct sensors to signal the 
appropriate type and amount of numerical-dissipation/filter for non-dissipative 
high order schemes. These scheme-independent sensors are capable of dis- 
tinguishing shocks/shears, turbulent fluctuations and spurious high-frequency 
oscillations. In addition, these sensors are readily available as an improvement 
over existing grid adaptation indicators. The same shock/shear detector that 
is designed to switch on the shock/shear numerical dissipation can be used to 
switch off the entropy splitting form of the inviscid flux derivative in the vicinity 
the discontinuous regions to further improve nonlinear stability and minimize 
the use of numerical dissipation. The rest of the sensors in conjunction with 
the local flow speed and Reynolds number can also be used to adaptively de- 
termine the appropriate entropy splitting parameter for each flow type/region. 
The minimization of employing very fine grids to overcome the production 
of spurious numerical solutions and/or instability due to under-resolved grids 
is also illustrated [79, 17]. Test examples shown are very encouraging. Full 
implementation of the approach for practical problems is in progress. 
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Appendix 

We here give boundary modified numerical dissipation operators, such that 
the semi bounded property is satisfied for the total operator. The operators 
are derived by splitting the periodic operator into one boundary part and one 
interior part, 

(D+D-) per = ( D + D_) b + (D+D-)i. (A.l) 

The interior operator (D + ) / has zeros in all rows where the centered operator 

can not be applied. We define the boundary modified operator as 


(D + D-y = (D+D-y^iD+D-f, 
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if q — 2p is even. The zero terms on the boundary of the interior operator, 
will destroy the periodic terms, so that the boundary modified operator does not 
have any periodic wrap around terms. If q = 2p + 1 is odd, we define 

{D+D.y = {D+D-f^D+^D-WD+D-Yj. 

The semi boundedness then follows by applying the summation by parts prop- 
erty for the periodic operator, or for the case of q even, 

(uj, ( D+DJ) q Uj) h = -(( D + D-)P er Uj , (D + D-) P jUj)h = 
-((D+D-YjUj^D+D-fjUjh^ 0 

where the last equality follows since the interior operator in the scalar product 
kills all boundary terms in the periodic operator. The case with odd q can be 
treated similarly. 

The order of the operator is reduced on the boundary. To have correct scaling, 
as numerical dissipation, the 2qlh derivative should be multiplied by h 2q ~ 1 , thus 
affecting the accuracy up to order 2q — 1. However, the boundary terms will 
be reduced to qth derivatives, and thus have order q- 1 on the boundary. For 
example, the fourth order operator will be third order in the interior, and first 
order on the boundary. The sixth order operator will be fifth order in the interior, 
and second order on the boundary, etc. We present below examples for orders 
4, 6, and 8. 
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Eighth derivative 
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